Graphics system

ABSTRACT

Shading values in a graphics processing system are determined by using approximation techniques. For example, where a Beckmann Distribution Function (BDF) is used, the shading values are determined by reformulating the BDF in terms of the scalar product of the surface normal and the viewpoint-light half-vector to produce a modified BDF. A portion of the BDF is divided into a plurality of intervals, each defined by first and second control points at its boundaries. Each interval is then approximated using a cubic approximation to produce coefficients for use in determining the shading values. In alternative embodiments, other functions are approximated.

[0001] The invention relates to data processing in a graphics system, and in particular, to real-time physically-based shading on SIMD (Single Instruction Multiple Data) architectures.

BACKGROUND OF THE INVENTION

[0002] With the advent of 3D graphics processors capable of per-pixel shading effects, focus has been shifting towards providing more advanced rendering techniques.

[0003] In recent years, much emphasis has been placed upon improving the speed of 3D graphics hardware. The focus has been to increase the number of polygons, the number of texels and so on. More recently, a trend has emerged to attempt to improve the quality of interactive 3D graphics, for example, by improving rendering techniques.

[0004] Traditionally, photo-realistic rendering has been the preserve of off-line rendering systems. However, more recently, these techniques have started to be used in contemporary interactive systems. Many of these techniques use the power of the texturing hardware to provide the desired functionality. For example, Cabral, et al., “Reflection Space Image Based Rendering, SIGGRAPH 99, conference proceedings, pages 165-170, and Heidrich et al., “Realistic, Hardware-accelerated Shading and Lighting, SIGGRAPH 99, Conference proceedings, pages 171-178, both demonstrate how effective this can be. The Cabral technique implements Bidirectional Reflection Distribution Functions (BRDF) using environment mapping hardware. Heidrich, on the other hand, uses multipass methods supported by current hardware to implement physically-based reflection models.

[0005] One example of a well known reflection model of a surface is the Phong model.

[0006] Another example of a well known physically-based reflection model is the Cook-Torrance model. This model is based on a previous model, known as the Torrance-Sparrow model, which describes a method that accounts for directional distribution of reflected light from roughened surfaces, simulated using microfacets.

[0007] Microfacets are assumed to be an isotropic collection of planar microscopic facets, each defined as a perfectly smooth reflector.

[0008] For the Cook-Torrance model, the specular component of bi-directional reflectivity is expressed as: $\rho_{s} = {\frac{F_{\lambda}}{\pi}\frac{DG}{\left( {\overset{\_}{N} \cdot \overset{\_}{V}} \right)\left( {\overset{\_}{N} \cdot \overset{\_}{L}} \right)}}$

[0009] where, F_(λ)=Fresnel term, relating incident light to reflected light

[0010] D=Distribution function for microfacets

[0011] G=Geometrical attenuation factor, the masking and shadowing effects of the microfacets on each other.

[0012] {overscore (N)}=Surface normal

[0013] {overscore (V)}=Viewpoint vector

[0014] {overscore (L)}=Light vector

[0015] In the above model, several known functions can be used to provide the microfacet distribution function, D. These include a Gaussian distribution; a Trowbridge-Reitz distribution and a Beckmann distribution. The Beckmann distribution is considered to be particularly advantageous because it provides the absolute magnitude of the reflectance without introducing arbitrary constants.

[0016] The Beckmann distribution function for rough surfaces is given as: $D = {\frac{1}{4m^{2}\cos^{4}\beta}^{- {(\frac{\tan \quad \beta}{m})}^{2}}}$

[0017] where, m=Root mean square slope of microfacets

[0018] β=Angle between {overscore (N)} and {overscore (H)}, where {overscore (H)} is the half vector between {overscore (V)} and {overscore (L)}.

[0019] The value of “m” is basically an indication of the smoothness or roughness of a surface. It ranges between 0 and 1, with smooth surfaces having low values of m, and rough surfaces having high values of m.

[0020] The Beckmann distribution function suffers from the disadvantage that it is not particularly suited for real-time implementation, since it is relatively expensive to compute, and it must be computed for each pixel that is being considered for shading. In particular, the cosine and tan functions require a relatively large number of instruction cycles to compute. This means that the Beckmann distribution function is not particularly suited for use on a SIMD architecture.

[0021] The computation power required to compute the Cook-Torrance reflectance model is further hampered by the Fresnel term, F, which is also expensive in terms of the number of instruction cycles required to calculate the term.

[0022] The Fresnel equation relates incident light to reflected light for the smooth surface of each microfacet. The Fresnel equation can be expressed as: $F_{\lambda} = {\frac{1}{2}\frac{1\left( {g - c} \right)^{2}}{2\left( {g + c} \right)^{2}}\left( {1 + \frac{\left( {{c\left( {g + c} \right)} - 1} \right)^{2}}{\left( {{c\left( {g - c} \right)} + 1} \right)^{2}}} \right)}$

[0023] where,

c=cosθ_(i) ·{overscore (L)}·{overscore (H)}

g ²=η_(λ) ² +c ²−1

[0024] $\eta_{\lambda} = \frac{\eta_{t}\lambda}{\eta_{i}\lambda}$

[0025] where, η_(t)λ and η_(i)λ are the indices of refraction of the two media.

[0026] As with the Beckmann distribution function mentioned above, the Fresnel term is also unsuited for real time calculations.

[0027] Thus, although known methods exist for offering improved shading techniques, they are not particularly suited for real-time calculations, or for use on a conventional SIMD architecture.

[0028] Therefore, the aim of the present invention is to provide an approximation technique, which can be used to provide a physically-based shading technique that requires less computation than known techniques, and which is suited for use on a SIMD architecture.

SUMMARY OF THE INVENTION

[0029] According to a first aspect of the present invention, there is provided method of determining shading values in a graphics processing system using an approximation of a Beckmann Distribution Function (BDF), the BDF being defined as: $D = {\frac{1}{4m^{2}\cos^{4}\beta}^{- {(\frac{\tan \quad \beta}{m})}^{2}}}$

[0030] where, m=Root mean square slope of microfacets

[0031] β=Angle between {overscore (N)} and {overscore (H)}

[0032] the method comprising the steps of:

[0033] (a) reformulating the BDF in terms of the scalar product of the surface normal {overscore (N)} and the viewpoint-light half-vector {overscore (H)} ({overscore (N)}·{overscore (H)}) to produce a modified BDF;

[0034] (b) dividing a portion of the BDF into a plurality of intervals, n, each interval defined by first and second control points at its boundaries; and,

[0035] (c) approximating each interval using a cubic approximation to produce coefficients for use in determining the shading values.

[0036] According to another aspect of the present invention, there is provided a method of determining shading values in a graphics processing system using a Cook-Torrance model, the Cook Torrance model being: $\rho_{s} = {\frac{F_{\lambda}}{\pi}\frac{DG}{\left( {\overset{\_}{N} \cdot \overset{\_}{V}} \right)\left( {\overset{\_}{N} \cdot \overset{\_}{L}} \right)}}$

[0037] where, F_(λ)=Fresnel term, relating incident light to reflected light

[0038] D=Distribution function for microfacets

[0039] G=Geometrical attenuation factor, the masking and shadowing effects of the microfacets on each other.

[0040] {overscore (N)}=Surface normal

[0041] {overscore (V)}=Viewpoint vector

[0042] {overscore (L)}=Light vector

[0043] the method comprising the steps of:

[0044] receiving data for a material being modelled;

[0045] approximating the microfacet distribution function D using a piecewise cubic interpolation of the distribution function; and,

[0046] calculating the Cook-Torrance model using the approximation of the microfacet distribution function.

[0047] According to another aspect of the present invention, there is provided a method of determining shading values in a graphics processing system using an approximation of a microfacet distribution function, the method comprising the steps of:

[0048] (a) dividing a portion of the microfacet distribution function into one or more intervals, n, each interval defined by first and second control points at its boundaries; and,

[0049] (b) approximating each interval using a polynomial approximation to produce coefficients for use in determining the shading values.

[0050] According to another aspect of the invention, there is provided a graphics processing system having means for determining shading values using an approximation of a Beckmann Distribution Function (BDF), the BDF being defined as: $D = {\frac{1}{4m^{2}\cos^{4}\beta}^{- {(\frac{\tan \quad \beta}{m})}^{2}}}$

[0051] where, m=Root mean square slope of microfacets

[0052] β=Angle between {overscore (N)} and {overscore (H)}

[0053] the graphics processing system comprising:

[0054] (a) reformulating means for reformulating the BDF in terms of the scalar product of the surface normal and the viewpoint-light half-vector ({overscore (N)}·{overscore (H)}) to produce a modified BDF;

[0055] (b) dividing means for dividing a portion of the BDF into a plurality of intervals, n, each interval defined by first and second control points at its boundaries; and,

[0056] (c) approximating means for approximating each interval using a cubic approximation to produce coefficients for use in determining the shading values.

[0057] According to another aspect of the invention, there is provided a graphics processing system for determining shading values using a Cook-Torrance model, the Cook Torrance model being: $\rho_{s} = {\frac{F_{\lambda}}{\pi}\frac{DG}{\left( {\overset{\_}{N} \cdot \overset{\_}{V}} \right)\left( {\overset{\_}{N} \cdot \overset{\_}{L}} \right)}}$

[0058] where, F_(λ)=Fresnel term, relating incident light to reflected light

[0059] D=Distribution function for microfacets

[0060] G=Geometrical attenuation factor, the masking and shadowing effects of the microfacets on each other.

[0061] {overscore (N)}=Surface normal

[0062] {overscore (V)}=Viewpoint vector

[0063] {overscore (L)}=Light vector

[0064] the graphics processing system comprising:

[0065] receiving means for receiving data for a material being modelled;

[0066] approximating means for approximating the microfacet distribution function D using piecewise cubic interpolation of the distribution function; and,

[0067] calculating means for calculating the Cook-Torrance model using the approximation of the microfacet distribution function.

[0068] According to another aspect of the present invention, there is provided a graphics processing system having means for determining shading values using an approximation of a microfacet distribution function, the graphics processing system comprising:

[0069] (a) dividing means for dividing a portion of the microfacet distribution function into one or more intervals, n, each interval defined by first and second control points at its boundaries; and,

[0070] (b) approximating means for approximating each interval using a polynomial approximation to produce coefficients for use in determining the shading values.

BRIEF DESCRIPTION OF THE DRAWINGS

[0071] For a better understanding of the present invention reference will now be made, by way of example, to the accompanying drawings, in which:

[0072]FIG. 1 shows a typical SIMD architecture;

[0073]FIG. 2 shows the steps involved in approximating a distribution function, such as the Beckmann distribution function, according to a first aspect of the present invention;

[0074]FIG. 3 shows how the approximation technique of FIG. 2 can be used in Cook-Torrance shading;

[0075]FIG. 4 shows Cook-Torrance shading having full Beckmann Distribution Function, and using the full Fresnel term, according to the prior art;

[0076]FIG. 5 shows Cook-Torrance shading having approximated Beckmann Distribution Function according to the present invention, plus an approximated Fresnel term.

DETAILED DESCRIPTION OF A PREFERRED EMBODIMENT OF THE INVENTION

[0077] The approximation techniques of the present invention will be described in relation to a SIMD based graphics hardware that can render scenes using physically-based reflection models. This is made possible by the ability of these systems to provide per-pixel calculation during the rendering process. However, it is noted that the approximation techniques are not limited for use of a SIMD architecture, but can be used in any architecture where the processing speed needs to be increased.

[0078]FIG. 1 shows a typical SIMD architecture, comprising a large array of general purpose ALU's, 1, connected using 1-D topology. The ALUs 1 are relatively simple, being 8-bit ALUs, for example, with limited floating point support. A similar architecture is described in the applicant's earlier PCT patent application no. PCT/GB00/01332.

[0079] Preferably, the SIMD architecture uses an OpenGL driver for medium to high level workstations that is used as the basis for adding Cook-Torrance shading as an extension. The processors 1 are preferably fully programmable, with 1536 8-bit processors arranged in 6 blocks of 256, 3₁ to 3₆. Preferably, a high bandwidth communication bus 5 exists to access the shared memory for the SIMD array.

[0080] The SIMD array has a graphics pipeline, provided by a software application, which is used for generating the graphics images. The lighting and shading are performed per pixel, resulting in unique surface normals, view vectors and light vectors.

[0081]FIG. 2 shows the steps involved in approximating a distribution function, such as the Beckmann distribution function, for use in such a SIMD architecture. It is noted that the technique is suitable for use with other distribution functions, such as Gaussian distribution and Trowbridge-Reitz distribution.

[0082] In order to reduce the amount of computation that is required each time the Beckmann distribution function is evaluated, the scalar product {overscore (N)}·{overscore (H)} is first calculated.

[0083] {overscore (N)} and {overscore (H)} are unit direction vectors which can be represented by cosβ. (H being a halfway vector between {overscore (V)} and {overscore (L)}).

[0084] The BDF is first reformulated in terms of {overscore (N)}·{overscore (H)}, step 1. This reformulated Beckmann Distribution Function will hereinafter be referred to a Beckmann_(mod), which is expressed as: $D = {\frac{1}{4{m^{2}\left( {\overset{\_}{N} \cdot \overset{\_}{H}} \right)}^{4}}^{{(\frac{\tan \quad \arccos \quad \beta}{m})}^{2}}}$

[0085] A portion of the Beckmann_(mod) curve is then divided into a plurality of intervals, n, as shown in step 3. Preferably, the portion of the curve relates to the most significant portion, i.e., the portion having the greatest change. The significant portion of the curve may be determined, for example, by detecting the portion of the curve where the value is non-zero. However, the portion of the curve could equally relate to the entire width of the curve.

[0086] Each of the n intervals will have endpoints, which are defined as the control points of that particular interval of the curve. Each interval is then approximated using a cubic equation, step 5. The use of a cubic equation has the advantage that smooth transitions are made possible near the intersection of two adjacent intervals. Discontinuity between adjacent intervals is avoiding by ensuring that the first order derivative is equal.

[0087] Preferably, the cubic equation is:

C ₁({overscore (N)}·{overscore (H)})³ +C ₂({overscore (N)}·{overscore (H)})² +C ₃({overscore (N)}·{overscore (H)})+C ₄

[0088] For each of the control points at either end of the interval in question, the coefficients C₁, C₂, C₃ and C₄ of the cubic equation are derived, step 7, using the following matrix: $\begin{bmatrix} C_{1} \\ C_{2} \\ C_{3} \\ C_{4} \end{bmatrix} = {\begin{bmatrix} x_{0}^{3} & x_{0}^{2} & x_{0} & 1 \\ x_{1}^{3} & x_{1}^{2} & x_{1} & 1 \\ {3x_{0}^{2}} & {2x_{0}} & 1 & 0 \\ {3x_{1}^{2}} & {2x_{1}} & 1 & 0 \end{bmatrix}^{- 1}\begin{bmatrix} P_{1} \\ P_{4} \\ R_{1} \\ R_{4} \end{bmatrix}}$

[0089] where x₀ and x₁ are the values of {overscore (N)}·{overscore (H)} for the interval endpoints; P₁ and P₄ are the values of Beckmann_(mod) at the endpoints and R₁ and R₄ are the values of the first order derivative at the interval endpoints (tangent to the function at x₀ and x₁).

[0090] The coefficients are then scaled by a factor of 2^(f), step 9, to bring them into the range −1 to 1. The value of “f” is chosen to be common to all control points. Therefore the max(abs(control point)) value is chosen, and f is then based on that value. The factoring of the coefficients in this manner allows the coefficients to be represented using 1.15 fixed point numbers, thereby avoiding the use of floating point calculations, which in turn provides a large saving on computation power (i.e. number of instructions required to compute the values).

[0091] A more detailed explanation will now be given regarding how the above mentioned matrix is derived.

[0092] The distribution function being approximated, for example the Beckmann distribution function, is relatively well behaved function in so far as its curve follows a relatively smooth path. As such, a piecewise approximation to the curve can give good results. The cubic piecewise approximation used in the present invention gives particularly good results, since it takes into account not only a series of control points, but also the tangent vectors at these points (i.e. the gradient of the curve, being the function's first order derivative).

[0093] Thus, the approximation technique calculates the first order derivative of the Beckmann function. Since the standard and first order derivatives of the function are being used in the interpolation of the curve, a Hermite curve can be used to carry out the approximation.

[0094] The following describes how the curve can be approximated in this way, using n control points. $\begin{matrix} {{{beckmann}\left( {\beta,m} \right)}:=\frac{e^{- {(\frac{\tan {(\beta)}}{m})}^{2}}}{\left( {4m^{2}{\cos (\beta)}} \right.}} \\ {{{dbeckmann}\left( {\beta,m} \right)}:={{- \left\lbrack {{2 \cdot \frac{\left\lbrack {1 + \left( {\tan (\beta)} \right)^{2}} \right\rbrack}{m^{2}{\cos (\beta)}}} + 1} \right\rbrack} \cdot e^{- {(\frac{\tan {(\beta)}}{m})}^{2}} \cdot \frac{\left. \left( {\tan \quad \beta} \right) \right)}{4m^{2}}}} \end{matrix}$ $\begin{matrix} {n:=\quad 8} \\ {{t:=\quad 0},{0.1\quad \ldots \quad n}} \\ {{{int}(x)}:=\quad {x - {{mod}\left( {x,1} \right)}}} \\ {{{low}(x)}:=\quad {{{int}(x)} \cdot \frac{u}{n}}} \\ {{{P1}(x)}:=\quad {{beckmann}\left( {{{low}(x)},m} \right)}} \\ {{{R1}(x)}:=\quad \frac{{debeckmann}\left( {{{low}(x)},m} \right)}{M}} \end{matrix}\begin{matrix} {\quad {u:=\frac{\pi}{2}}} \\ \quad \\ \quad \end{matrix}\begin{matrix} {m:=\quad 0.2} \\ {M:=\quad 4} \\ \quad \\ {{{hi}(x)}:=\quad \left( {{{int}(x)} + {1 \cdot \frac{u}{n}}} \right.} \\ {{{P4}(x)}:=\quad {{beckmann}\left( {{{hi}(x)},m} \right)}} \\ {{{R4}(x)}:=\quad \frac{{dbeckmann}\left( {{{hi}(x)},m} \right)}{M}} \end{matrix}$ $\begin{matrix} {{Q\left( {u,{p1},{p4},{r1},{r4}} \right)}:=\quad {{\left( {{2u^{3}} - {3u^{2}} + 1} \right){p1}} + {\left( {{{- 2}u^{3}} + {3u^{2}}} \right){p4}} +}} \\ {\quad {{\left( {u^{3} - {2u^{2}} + u} \right){r1}} + {\left( {u^{3} - u^{2}} \right){r4}}}} \\ {{{hermy}(u)}:=\quad {Q\left( {{u - {{int}(u)}},{{P1}(u)},{{P4}(u)},{{R1}(u)},{{R4}(u)}} \right)}} \\ {{{hermx}(u)}:=\quad {Q\left( {{u - {{int}(u)}},{{low}(u)},{{hi}(u)},\frac{1}{M},\frac{1}{M}} \right)}} \end{matrix}$

[0095] The Hermite curves are parametric curves defined in t. Thus, when plotting the Beckmann curve as an approximation using Hermite curves, it is necessary to plot two functions on the axes of the graph. In the example described above, it is necessary to plot hermy(t) against hermx(t). This has the advantage that for a given x there can be multiple values of y. In essence, this means that loops can be created in the curve. For well behaved curves, this is unnecessary since the curve should not contain loops. For well behaved curves, the function defining x should be very close to a multiple of t. However, in the case of the Beckmann function, the curve is not well behaved when replacing kt for hermx(t). It is approximately the same curve, but is not as smooth. This results in mach bands appearing on the eventual image. In order to make the curve appear smooth, it is necessary to use a different class of curve, which is defined around x directly instead of t. This is explained more fully in the following set of equations.

[0096] The hermite curve is defined as the following cubic equation defined using parametric term t.

Q(t)=C1t ³ +C2t ² +C3t+C4 where 0.0≦t≦1.0

[0097] There are defined a set of coefficients (C₁, C₂, C₃ and C₄) defined for each interval. The endpoints of the interval are defined as the control points of the curve.

[0098] Furthermore, we can define the function as follows

x(t)=C1_(x) t ³ +C2_(x) t ² +C3_(x) t+C4_(x)

[0099] Where C_(1x), C_(2x), C_(3x), C_(4x) are the x components of the coefficients.

[0100] By differentiating the equation shown above, and substituting values for t at t=0 and t=1 (where the known control points exist), it is possible to create a matrix, representing the equations for the control points, and their derivatives such that $\begin{bmatrix} P_{1} \\ P_{4} \\ R_{1} \\ R_{4} \end{bmatrix} = {{G_{Hx}\left\lbrack \quad \begin{matrix} 0 & 0 & 0 & 1 \\ 1 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 \\ 3 & 2 & 1 & 0 \end{matrix}\quad \right\rbrack} \cdot M_{H} \cdot G_{Hx}}$

[0101] where the rows of the 4×4 matrix are simply Q(t) with t=0 and t=1 substituted (rows 0 and 1), and Q′(t) with t=0 and t=1 substituted (rows 2 and 3). M_(H) is the matrix which is used to multiply the control points to give the coefficients. P₁, P₄, R₁ and R₄ are the control points at t=0 and t=1 and the tangent vectors to the curve at t=0 and t=1 (derived in this case by differentiating the Beckmann function and substituting in the relevant values).

[0102] Thus we can derive a value for M_(H) as follows $M_{H} = {\left\lbrack \quad \begin{matrix} 0 & 0 & 0 & 1 \\ 1 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 \\ 3 & 2 & 1 & 0 \end{matrix}\quad \right\rbrack^{- 1} = \left\lbrack \quad \begin{matrix} 2 & {- 2} & 1 & 1 \\ {- 3} & 3 & {- 2} & 1 \\ 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 \end{matrix}\quad \right\rbrack}$

[0103] which allows us to define the parametric form of the Hermite curve in terms of the control points as follows

Q(t)+(2t ³−3t ²+1)P ₁+(−2t ³+3t ²)P ₄+(t ³−2t ² +t)R ₁+(t ³ −t ²)R ₄

[0104] It is clear from this curve that for values of t=0, the value of Q (t) is P₁ which is exactly what one would expect.

[0105] As indicated earlier, this form of the curve requires the function to be evaluated for both x and y. Using an approximation of x based on multiples of t is a reasonable approximation, but does mean that mach bands appear where the curve is not as smooth as the fully derived case. In order to combat this, the curve is redefined to be derived as follows.

Q(x)=C1x ³ +C2x ² +C3x+C4

[0106] This allows a smooth curve to be defined for the function for a given x. This introduces several problems. Since the interval between control points is no longer defined over 0 to 1, a standard matrix to derive the coefficients cannot be used. Since the curve is defined over x, it is not possible to have multiple values of y for a given x. This is not a problem for well behaved curves, but may pose a problem for curves that are not well behaved, since the curve will be undefined.

[0107] In order to define the new function, it is necessary to reformulate the matrix and its inverse. Using the same values as control points on the curve as before, i.e P₁ and P₄ (points on the curve) and R₁ and R₄ (tangents to the curve at the control points), the coefficient definition can be reformulated as follows. $\begin{bmatrix} C_{1} \\ C_{2} \\ C_{3} \\ C_{4} \end{bmatrix} = {\left\lbrack \quad \begin{matrix} x_{0}^{3} & x_{0}^{2} & x_{o} & 1 \\ x_{1}^{3} & x_{1}^{2} & x_{1} & 1 \\ {3x_{o}^{2}} & {2x_{0}} & 1 & 0 \\ {3x_{1}^{2}} & {2x_{1}} & 1 & 0 \end{matrix}\quad \right\rbrack^{- 1}\begin{bmatrix} P_{1} \\ P_{4} \\ R_{1} \\ R_{4} \end{bmatrix}}$

[0108] where x₀ and x₁ are the x coordinates at the endpoints of the interval.

[0109] Thus, the matrix must be created and inverted for each curve interval, which gives the four coefficients for the cubic curve equation.

[0110] This produces a well behaved curve which is defined for x.

[0111] The above mentioned approximation provides a very close match to the true Beckmann distribution function, but requires far less computation. In addition, by factoring the results by a factor of 2^(f), the coefficients can be normalised to a range −1.0 to 1.0, thereby avoiding the need for floating point calculations. Instead, the final result can be achieved by shifting the result by the power of the factor value, f. This shift technique, which can either be a shift-left or a shift-right requires minimal computational work.

[0112] Each modified Beckmann distribution function, Beckmann_(mod), can be represented by a 56 byte look-up table. This table can either be embedded in the material for an object, or can be calculated when a material is loaded into the context for a scene, i.e. calculated on the fly.

[0113] For this purpose, each processor element in the SIMD array must be able to access a memory location using an offset unique to that element. This enables the array to access look-up tables in parallel across the array.

[0114] Once the coefficients have been extracted from the look up table, the cubic equation is multiplied out, for example, by expressing as 3 multiplies and 3 additions. This is a significant saving over the original Beckmann Distribution Function calculation.

[0115] The approximation technique mentioned above can be applied to distribution functions other than the Beckmann function. In addition to other microfacet distribution functions such as the Gaussian and Trowbridge-Reitz distribution functions, the approximation technique could also be used to simplify the calculation of the Fresnel term.

[0116] A piecewise cubic approximation is formed for the Fresnel curve for given wavelengths, and a suitable look-up table created in the material, in a similar manner to that already described above in relation to the Beckmann function.

[0117] Alternatively, when calculating the Cook-Torrance model, the Fresnel term could be approximated using other techniques. For example, Cook and Torrance use an approximation for each of the Red, Green and Blue components. According to another example, the calculation of Cook-Torrance model could use a simpler approximation for the Fresnel term, which assumes that the viewpoint and light source are at infinity. This simplifies matters in that the equation need only be calculated once per frame, as it will be the same at all points in the scene. However, this has a disadvantage is that inaccuracies can become apparent as the light and/or viewer move close to a surface.

[0118] Referring to FIG. 3, the Cook-Torrance model calculates the shading based on the approximations of the Beckmann function and/or Fresnel term mentioned above.

[0119] The SIMD architecture receives data relating to the material which is being modelled, step 1. During the lighting stage in the graphics pipeline, at each pixel, various other values in the Cook-Torrance equation will have been calculated. Thus, most of the terms may have already been calculated in other parts of the lighting calculation, and are therefore already available for the Cook-Torrance model, step 3. These include, for example, the surface normal N, the viewpoint vector V, the light vector L, and H, the halfway vector between V and L. The value of m (i.e. the root mean square slope of microfacets) will be embedded in the material, and will be constant for a given material.

[0120] The Beckmann_(mod) data may either be available in the form of a look-up table embedded in the material, or, if not, is calculated once per material. Once calculated, the Beckmann_(mod) data is then available for subsequent data processing, and used to approximate the microfacet distribution function, step 5.

[0121] If an approximation of the Fresnel term is also being used, then this approximation is derived, steps 7, 9. Otherwise, the full Fresnel tern is calculated, step 11.

[0122] The Cook-Torrance model is then calculated, step 13, using the approximation(s) previously calculated.

[0123] It is noted that although the preferred embodiment refers to using the Beckmann Distribution Function, the invention is equally applicable to any other microfacet distribution function. Likewise, although the preferred approximation of the Fresnel term uses a piecewise cubic approximation similar to the Beckmann distribution function, other approximations of the Fresnel term may also be used. Furthermore, the Fresnel term may be approximated using an array wide value which has already been calculated.

[0124] Also, although the preferred embodiment uses a cubic approximation, other polynomials may be used, such a linear, quadratic, quintic, and so on.

[0125] The approximation methods described above enable the Cook-Torrance model to be computed quicker, thereby allowing real-time physically-based shading to be performed. The method of data processing is particularly suited for use in a SIMD environment.

[0126] FIGS. 4 illustrates a set of spheres rendered using a Cook-Torrance model using the full Beckmann distribution function and full Fresnel term, i.e. no approximation. The top eight spheres use m=0.2, while the bottom eight spheres use m=0.6. FIG. 5 shows the spheres rendered using a Cook-Torrance shading using the piecewise cubic approximation of the present invention, again having the top eight spheres using m=0.2 and the bottom eight spheres using m=0.6, but where the Fresnel term was approximated using infinite light source and viewpoint. This approximation of the Fresnel term accounts for the slight differences in the results. The light source and viewpoint in these images were quite close to the surface of the sphere in order to accentuate any problems this approximation would cause. In most cases, the light source and viewpoint would be much further from the object and thus cause less marked differences. However, even in this extreme case, the difference is not significant. Without any comparison with the correct result an observer would not notice the small error.

[0127] Using the approximations described above, the Cook-Torrance model can be implemented on a SIMD architecture in approximately 1000 cycles. If, for example, the architecture uses 1536 processors running at 200 MHZ, this translates to 307 Mpixels/s. This compares very favourably with a similar architecture implementing Phong shading, which also takes approximately 1000 cycles.

[0128] If we take into account the fact that the hardware must also perform all geometry calculation, texturing, diffuse lighting calculations and other features in the rendering pipeline for OpenGL, it is evident that the Phong and Cook-Torrance terms are a minor part of the overall time taken per frame. This reduces the effective cost of the Cook-Torrance shading relative to the conventional Phong shading. 

What is claimed is:
 1. A method of determining shading values in a graphics processing system using an approximation of a Beckmann Distribution Function (BDF), the BDF being defined as: $D = {\frac{1}{4m^{2}\cos^{4}\beta}e^{- {(\frac{\tan \quad \beta}{m})}^{2}}}$

where, m=Root mean square slope of microfacets β=Angle between {overscore (N)} and {overscore (H)} the method comprising the steps of: (a) reformulating the BDF in terms of the scalar product of the surface normal {overscore (N)} and the viewpoint-light half-vector {overscore (H)} ({overscore (N)}·{overscore (H)}) to produce a modified BDF; (b) dividing a portion of the BDF into a plurality of intervals, n, each interval defined by first and second control points at its boundaries; and, (c) approximating each interval using a cubic approximation to produce coefficients for use in determining the shading values.
 2. A method as claimed in claim 1, wherein the coefficients are scaled by a factor of 2^(f).
 3. A method as claimed in claim 2, whereby f is equal for all intervals.
 4. A method as claimed in claim 1, wherein the modified Beckmann distribution function is defined as: $D = {\frac{1}{4\left( {m^{2}\left( {\overset{\_}{N} \cdot \overset{\_}{H}} \right)} \right)^{4}}{e^{- {(\frac{\tan \quad {arc}\quad \cos \quad \beta}{m})}^{2}}.}}$


5. A method as claimed in claim 1, wherein the step of approximating each interval using a cubic approximation uses the following matrix to determine the coefficients of the cubic equation: $\begin{bmatrix} C_{1} \\ C_{2} \\ C_{3} \\ C_{4} \end{bmatrix} = {\left\lbrack \quad \begin{matrix} x_{0}^{3} & x_{0}^{2} & x_{0} & 1 \\ x_{1}^{3} & x_{1}^{2} & x_{1} & 1 \\ {3x_{0}^{2}} & {2x_{0}} & 1 & 0 \\ {3x_{1}^{2}} & {2x_{1}} & 1 & 0 \end{matrix}\quad \right\rbrack^{- 1}\begin{bmatrix} P_{1} \\ P_{4} \\ R_{1} \\ R_{4} \end{bmatrix}}$

where, P=value of BDF at the end points of the interval, and R=value of 1^(st) order derivatives at the end points of the interval.
 6. A method as claimed in claim 1, wherein the portion of the BDF relates to a portion in which the data changes significantly.
 7. A method as claimed in claim 1, wherein the modified BDF data is calculated once for a given material.
 8. A method as claimed in claim 1, wherein the modified BDF data is predetermined for a given material.
 9. A method as claimed in claim 1, wherein the cubic equation is: C ₁({overscore (N)}·{overscore (H)})³ +C ₂({overscore (N)}·{overscore (H)})² +C ₃({overscore (N)}·{overscore (H)})+C ₄
 10. A method as claimed in claim 1, wherein n=8.
 11. A method of determining shading values in a graphics processing system using a Cook-Torrance model, the Cook-Torrance model being: $\rho_{s} = {\frac{F_{\lambda}}{\pi}\frac{DG}{\left( {\overset{\_}{N} \cdot \overset{\_}{V}} \right)\left( {\overset{\_}{N} \cdot \overset{\_}{L}} \right)}}$

where, F_(λ)=Fresnel term, relating incident light to reflected light D=Distribution function for microfacets G=Geometrical attenuation factor, the masking and shadowing effects of the microfacets on each other. {overscore (N)}=Surface normal {overscore (V)}=Viewpoint vector {overscore (L)}=Light vector the method comprising: receiving data for a material being modelled; approximating the microfacet distribution function D using a piecewise cubic interpolation of the distribution function; and, calculating the Cook-Torrance model using the approximation of the microfacet distribution function.
 12. A method as claimed in claim 11, wherein the step of approximating the microfacet distribution function comprises the steps of: determining whether a look-up table exists for the set of coefficients for approximating the BDF using cubic interpolation, and, if so, extracting the appropriate data, and, if not, calculating the data by: dividing the function into a plurality of intervals, n; approximating each interval using a cubic approximation; for each control point at the end of an interval, deriving the coefficients of the cubic equation; scaling the coefficients by a factor of 2^(f); and, multiplying out the cubic equation to obtain the microfacet distribution values.
 13. A method of determining shading values in a graphics processing system using an approximation of a microfacet distribution function, the method comprising: (a) dividing a portion of the microfacet distribution function into one or more intervals, n, each interval defined by first and second control points at its boundaries; and, (b) approximating each interval using a polynomial approximation to produce coefficients for use in determining the shading values.
 14. A method as claimed in claim 13, wherein the polynomial approximation is a linear approximation.
 15. A method as claimed in claim 13, wherein the polynomial approximation is a quadratic approximation.
 16. A method as claimed in claim 13, wherein the polynomial approximation is a cubic approximation.
 17. A method as claimed in claim 13, wherein the coefficients are scaled by a factor of 2^(f).
 18. A method as claimed in claim 17, wherein f is equal for all intervals.
 19. A method as claimed in claim 13, wherein the portion of the microfacet distribution function relates to a portion in which the data changes significantly.
 20. A method as claimed in claim 13, wherein n=8.
 21. A method as claimed in claim 13, wherein the shading values are determined in real time.
 22. A graphics processing system having means for determining shading values using an approximation of a Beckmann Distribution Function (BDF), the BDF being defined as: $D = {\frac{1}{4m^{2}\cos^{4}\beta}e^{- {(\frac{\tan \quad \beta}{m})}^{2}}}$

where, m=Root mean square slope of microfacets β=Angle between {overscore (N)} and {overscore (H)} the graphics processing system comprising: (a) reformulating means for reformulating the BDF in terms of the scalar product of the surface normal and the viewpoint-light half-vector ({overscore (N)}·{overscore (H)}) to produce a modified BDF; (b) dividing means for dividing a portion of the BDF into a plurality of intervals, n, each interval defined by first and second control points at its boundaries; and, (c) approximating means for approximating each interval using a cubic approximation to produce coefficients for use in determining the shading values.
 23. A graphics processing system as claimed in claim 22, further comprising factoring means for scaling the coefficients by a factor of 2^(f).
 24. A graphics processing system as claimed in claim 23, whereby f is equal for all intervals.
 25. A graphics processing system as claimed in claim 22, wherein the modified Beckmann distribution function is defined as: $D = {\frac{1}{4\left( {m^{2}\left( {\overset{\_}{N} \cdot \overset{\_}{H}} \right)} \right)^{4}}{e^{- {(\frac{\tan \quad {arc}\quad \cos \quad \beta}{m})}^{2}}.}}$


26. A graphics processing system as claimed in claim 22, wherein the approximating means comprises the following matrix to determine the coefficients of the cubic equation: $\begin{bmatrix} C_{1} \\ C_{2} \\ C_{3} \\ C_{4} \end{bmatrix} = {\begin{bmatrix} x_{0}^{3} & x_{0}^{2} & x_{0} & 1 \\ x_{1}^{3} & x_{1}^{2} & x_{1} & 1 \\ {3x_{0}^{2}} & {2x_{0}} & 1 & 0 \\ {3x_{1}^{2}} & {2x_{1}} & 1 & 0 \end{bmatrix}^{- 1}\begin{bmatrix} P_{1} \\ P_{4} \\ R_{1} \\ R_{4} \end{bmatrix}}$

where, P=value of BDF at the end points of the interval, and R=value of 1^(st) order derivatives at the end points of the interval.
 27. A graphics processing system as claimed in claim 22, wherein the dividing means divides a portion of the BDF in which the data changes significantly.
 28. A graphics processing system as claimed in claim 22, wherein the reformulating means calculates the modified BDF data once for a given material.
 29. A graphics processing system as claimed in claim 22, wherein the modified BDF data is predetermined for a given material.
 30. A graphics processing system as claimed in claim 22, wherein the approximating means uses the following cubic equation: C ₁({overscore (N)}·{overscore (H)})³ +C ₂({overscore (N)}·{overscore (H)})² +C ₃({overscore (N)}·{overscore (H)})+C ₄
 31. A graphics processing system as claimed in claim 22, wherein n=8.
 32. A graphics processing system for determining shading values using a Cook-Torrance model, the Cook Torrance model being: $\rho_{s} = {\frac{F_{\lambda}}{\pi}\frac{DG}{\left( {\overset{\_}{N} \cdot \overset{\_}{V}} \right)\left( {\overset{\_}{N} \cdot \overset{\_}{L}} \right)}}$

where, F_(λ)=Fresnel term, relating incident light to reflected light D=Distribution function for microfacets G=Geometrical attenuation factor, the masking and shadowing effects of the microfacets on each other. {overscore (N)}=Surface normal {overscore (V)}=Viewpoint vector {overscore (L)}=Light vector the graphics processing system comprising: receiving means for receiving data for a material being modelled; approximating means for approximating the microfacet distribution function D using piecewise cubic interpolation of the distribution function; and, calculating means for calculating the Cook-Torrance model using the approximation of the microfacet distribution function.
 33. A graphics processing system as claimed in claim 32, wherein the approximating means comprises: means for determining whether a look-up table exists for the set of coefficients for approximating the microfacet distribution function using cubic interpolation, and, if so, extracting the appropriate data, and, if not, calculating the data by: dividing the function into a plurality of intervals, n; approximating each interval using a cubic approximation; for each control point at the end of an interval, deriving the coefficients of the cubic equation; scaling the coefficients by a factor of 2^(f); and, multiplying out the cubic equation to obtain the microfacet distribution value.
 34. A graphics processing system having means for determining shading values using an approximation of a microfacet distribution function, the graphics processing system comprising: (a) dividing means for dividing a portion of the microfacet distribution function into one or more intervals, n, each interval defined by first and second control points at its boundaries; and, (b) approximating means for approximating each interval using a polynomial approximation to produce coefficients for use in determining the shading values.
 35. A graphics processing system as claimed in claim 34, wherein the approximating means uses a linear approximation.
 36. A graphics processing system as claimed in claim 34, wherein the approximating means uses a quadratic approximation.
 37. A graphics processing system as claimed in claim 34, wherein the approximating means uses a cubic approximation.
 38. A graphics processing system as claimed in claim 32, wherein the coefficients are scaled by a factor of 2^(f).
 39. A graphics processing system as claimed in claim 38, wherein f is equal for all intervals.
 40. A graphics processing system as claimed in claim 34, wherein the portion of the microfacet distribution function relates to a portion in which the data changes significantly.
 41. A graphics processing system as claimed in claim 32, wherein n=8.
 42. A graphics processing system as claimed in claim 34, wherein the coefficients are scaled by a factor of 2^(f).
 43. A graphics processing system as claimed in claim 34, wherein n=8. 